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Abstract 

Recent models for discrete euclidean quantum gravity incorporate a sum over 
simplicial triangulations. We describe an algorithm for simulating such models in 
arbitrary dimension. As illustration we show results from simulations in four dimen- 
sions. 
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Introduction 



In recent years there has been considerable interest in studying statistical systems whose 
partition functions encorporate a sum over simplicial triangulations. Initially, efforts fo- 
cused on two dimensional models which were proposed as discrete regularisations for string 
theory out of the critical dimension |], H, §]. Plausibility arguments were given to suggest 
that the sum over lattices, in some scaling limit, would generate the effects of a nonper- 
turbative inclusion of fluctuations in the worldsheet geometry Later, a framework was 
developed 0, |5| which, in certain simple cases, allowed many of the important features of 
a system coupled to two dimensional gravity to be derived from a knowledge of the theory 
in flat space. 

The results of these continuum analyses were in complete agreement with calculations 
and numerical simulations of the triangulated models || [?|, §, |9], [H], |TTJ and lent strong 
support to the lattice prescription. The continuum methods used to analyse the simple 
models (including pure two dimensional gravity) appear to break down for strings in 
physical dimensions. This has motivated a variety of numerical studies of the discrete 



models (which are well defined everywhere) with some interesting results [12|, [13|, [14|, [15 

In addition, the basic idea of summing over simplicial triangulations to generate a 
path integral for quantum gravity has been extended to three 0, |17], [Tj| [H| and four 



dimensions [EUl BTL E3L pi . Whilst these initial simulations are rather exploratory, the 



results for four dimensions are particularly exciting, as they seem to hint at a nonpertur- 
bative fixed point in the quantum theory. It is possible that the problems associated with 
the nonrenormalisability of Einstein gravity might be evaded in any continuum theory 
constructed in the vicinity of this new fixed point. 

In this paper we present an algorithm for Monte Carlo simulation of these dynamically 
triangulated models, which is constructed in such a way as to make trivial the dependence 
of the code on the manifold dimension d. 

Model 

We will be considering the problem of estimating a partition function of the form 

Z{Ko,K d )= e~ s(K0 ' K - T) (1) 

T(s d ) 

The summation goes over all simplicial triangulations T of the sphere in d dimensions 
S d . A simplicial triangulation is specified by a set of <i-simplices (sets of d + 1 labelled 
points) which are associated uniquely in pairs via their d — 1-dimensional faces, in such 
a way that the neighbourhood of any point is homeomorphic to a d- dimensional ball. In 
two dimensions d = 2, the fundamental building blocks are 2-simplices (triangles) which 
are glued together along their 1-dimensional faces (links) in such a way that two points 
(vertices) are connected by at most one link and the end points of all links are distinct. 
Similarly, a three dimensional simplicial manifold is built out of 3-simplices (tetrahedra) 
such that a given 2-dimensional face (triangle) is associated with exactly two tetrahedra. 
The restriction to manifolds ensures that every subsimplex is nondegenerate and unique. 
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Analogously, in four dimensions, the triangulation consists of hypertetrahedra associated 
in pairs via their tetrahedral faces. Again, the manifold condition effectively eliminates 
any degeneracies in the subsimplices. 

The action S for dimensions d < 4 can be taken to depend on only two coupling 
constants k , n d related to the bare Newton and cosmological constants respectively. They 
are conjugate to the total number of (i-simplices N d and points N (O-simplices) for a given 
triangulation T. 

S (/t , K d ) = K d N d - KqNq (2) 

The numerical evaluation of this partition function Z (and expectation values computed 
from it) is effected by a Monte Carlo procedure which generates a random walk in the 
space of all such triangulations by a sequence of local 'moves' or deformations. The ones 
commonly used correspond to the replacement of an z-dimensional subsimplex (i.e a subset 
of % + 1 points within a simplex) by its 'dual' d — i subsimplex (see, for example, [p5||). In 
order that this move preserve the manifold structure of the triangulation, there will be an 
associated change in the number and identity of neighbouring simplices and subsimplices. 
These moves have been shown to be ergodic (at least when d < 4) in ||25| . The latter 
statement implies that, at least in principle, a set of such moves can transform any such 
triangulation into any other of the same Euler character \. The topological invariant \ is 
defined for triangulations as 

X = E(~im (3) 

i=0 

However, in the case of d — 4, it appears that the typical number of such moves may 
increase very rapidly with volume (number of rf-simplices). This may place important 
constraints on the 'practical' ergodicity of the numerical simulations pBJ. Indeed, there 
is some recent numerical evidence that in this case the triangulation space may grow 
factorially with volume | ^7|| . 

Clearly, there are d + 1 possible moves which may be labeled by the dimension of the 
subsimplex i central to the move. The definition of the i-move requires that the order of 
a given z-subsimplex (the number of <i-simplices associated with it) be exactly d + 1 — i. 
Such a subsimplex will be referred to as a legal subsimplex. Subsequent upon finding such 
a legal subsimplex it is necessary to test whether substituting it by its 'dual' will lead to 
a bona fide triangulation which satisfies the manifold restriction. Such a move is referred 
to as geometrically allowed. Finally, for such a geometrically allowed move, the change in 
the action is computed and the update subjected to a metropolis test. Together with an 
explicit detailed balance condition, this procedure, under repeated iteration, will guarantee 
that the configurations approach a static distribution governed by the Boltzmann weight. 
The details of this algorithm are given in the next section. 

Algorithm 

The recipe for generating legal moves is as follows. Select a move type i at random. Then 
choose a simplex and one of its i-subsimplices ((i + 1) vertex labels) also at random. Using 



2 



a local procedure find the order of this subsimplex (i.e the number of simplices to which 
it is associated). If O (i) (d + 1 — i) go back and select another move type. 

The details of the neighbour search are organised as follows. Denote the labels of the 
(d + 1) points making up a simplex containing the subsimplex % in question by a . . . ad- 
Examine all neighbour simplices which are associated with this simplex by any face con- 
taining the i-subsimplex. There are d — i of these and each contains one vertex which is 
not in the original simplex. If the move is to be 'legal' (i.e the subsimplex i is of order 
d + 1 — i) then this extra vertex must be the same in all these d — i cases and can be de- 
noted ad+i. The slight exception to this picture corresponds to barycentric node insertion 
(i = d) where the extra vertex is a new label and no searching is required. 

It is now convenient to relabel the d + 2 vertices central to the move in such a way that 
the i + 1 points that define the subsimplex are arranged from to a iy the content of the 
i-move may then be seen from the following construction 

(4) 

The d + 1 — i initial state simplices (the lefthand side of this equation) are constructed 
by pairing the common subsimplex vertices (indicated by the brace) with d — i selected 
from the d — i + 1 other vertices. The final state is now gotten by identifying the points 
di + i . . . dddd+i as the new common subsimplex vertices (as indicated by the shift of the 
brace). The vertices needed to make up the % + 1 final state simplices are just % selected 
from the i + 1 remaining oq • • • a$ ■ 

In order that the new simplicial complex still corresponds to a triangulation of a mani- 
fold, it is necessary to check that the potential new simplices and subsimplices introduced 
by such a move are not already present in the triangulation. In effect, this means that the 
extra vertex dd+i must not already exist in any simplex associated to the subset dj+i . . . ad- 
To check for this a local search is carried out on all simplices which contain this subset. 
The nearby simplices are explored by moving out on faces containing this subset, with 
simplices being flagged and removed from the search list when they have been examined 
once. 

Once this manifold condition has been checked, the update is treated by the usual 
Metropolis test and the triangulation updated if necessary. In order that the simulation 
produce the correct Boltzmann probability density we have chosen to encorporate strict 
detailed balance into the algorithm (see |[23|| ). Denoting the probability of transition 
between one state or triangulation a to another (3 via some subsimplex move % by r (z, a, j3), 
detailed balance requires 

P(a)r(i,a,p) = P(p)r(d-i,P,a) (5) 

P (a) is the usual factor 

P ( a ) = e- si - a) (6) 

In practice, the transition probability factors into a product of probabilities to select the 
initial and final states - r\ (i, a) (f> (d — i, (3) , together with a piece t (i, a, (3) dependent on 
the change in action S. Here, a choice of initial state produces a unique final state so 
0=1. 
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In practice, an attempted update starts with a random selection of the move type i 
followed by a random selection of a simplex. The probability then of selecting a given 
i-subsimplex is then where O (i) is the order of the subsimplex in the triangulation 
a. For a legal move O (i) = {d + 1 — i). If the % + 1 vertices are then drawn at random 
from this simplex, the total probability of selection is 

V[h } d + lN d (a) [7) 

It is elementary to then see that the inverse move rj (d — i, j3) differs only by the number 
of simplices N d {(5). Thus eqn. 5 reads 

e ~ S{a) WT~^ & a > & = e ~ m WToS* V-t'M ( 8 ) 

This relation is then satisfied by choosing the reduced transition matrix t (i, a, (3) to have 
the simple form 

t (i, a, 3) = t j- (9) 

1+1+ ^~ d }) e SW-S(a) KJ 

The change in action only depends on the order % of the subsimplex move (since the total 
change in the number of simplices is 2i — d). 

S(J3)-S (a) = K d (2i -d)-n (5 iid - Sijo) +^(2t- d) (2 (N d — V) + 2i — d) (10) 

The supplementary term with coefficient 7 acts to control the volume fluctuations so that 
with a tuning of K d we can simulate a quasi microcanonical ensemble of fixed volume V. 
Typically the coupling 7 is taken small 7 = 0.005. We have verified that expectation 
values computed at small 7 are independent of 7 but possess statistical errors that grow 
as 7 — > 0. 

To simulate a lattice with volume V we tune the bare cosmological coupling n d during 
equilibration according to a formula which follows from steepest descent evaluation of the 
partition function 

5K d = 2 1 ((N d )-V) (11) 



Data structures and practical considerations 

The code is written in C in order to handle dynamic memory allocation. A structure 
of type SIMPLEX is defined which contains both an array of labels for its vertices and 
an array of pointers to the neighbour simplices. The pointer to a neighbour simplex is 
stored with a local array index identical to the vertex which does not appear in the face 
separating the two simplices. In addition each simplex contains a logical flag which may 
be set and unset during simplex searching operations to prevent a given simplex being 
used more than once. Finally the sum of its labels is also stored, as this allows for a fast 



calculation of the opposing vertex of a new simplex neighbour across a given face p3 . 
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Sequences of functions allocate and delete simplices dynamically and update the pointer 
fields of simplices neighbour to a move. To handle the node insertion and deletion moves, 
a stack of 'used' vertex labels is maintained. If a node insertion is attempted, the new 
label is drawn from the top of this stack, unless the stack is empty in which case the total 
node number is incremented. Conversely deleted nodes are placed on the stack. The stack 
itself is managed as a linked list. 

The total storage is of order 4 (2d + 12) V bytes for a V simplex simulation. This 
equates to approximately 0.6 Mbyte for an 8000 simplex lattice in four dimensions. 

The update time for V = 8000, d = 4 at Ko = 0.0 is of order 4000 microsecs per 
accepted elementary move on a HP-735 workstation. One Monte Carlo sweep is defined 
as V attempted, legal subsimplex moves. At V = 8000 our sweep time is 1.4 sees - 
this equates to 180 microsecs per attempted move and hence our average acceptance 
rate is approximately 5%. Sweep times decrease monotonically with increasing k due to 
the decreased connectivity of the lattice yielding significantly faster local search times. 
At k = 2.4 (close to criticality for V = 8000) the sweep time is just 0.9 sees with a 
correspondingly smaller update time per accepted move. 

The CPU time per attempted move increases with volume in d — 4 and for small Kq: 
for V = 4000 it is 120 microsec, at V = 8000 it is 180 microsec and for V = 16000 it has 
reached 230 microsec. However for sufficiently large kq it is essentially constant. 

Fig. [l] illustrates a typical execution profile for the code with d = 4 and V = 8000, 
giving the percentage CPU time spent in the most important routines which are labelled 
according to 

1. Searching for legal subsimplices 

2. Checking the geometric constraints 

3. Computing the metropolis test 

4. Updating the lattice structures 

Clearly, the program is dominated by the searching required to check that a proposed 
move does not violate the geometric restriction to manifolds. 

Characteristic output 

Fig. is a histogram illustrating A (i) the number of accepted moves of type i per sweep 
for a four dimensional lattice of volume V = 8000 at zero node coupling. The manifest 
symmetry about % = d/2 is a crude check of detailed balance - there are as many moves 
of type i as inverse moves of type d — i. 

For the same run, fig. [3] shows L (i) the average number of legal subsimplices of type 
i encountered per sweep. Clearly, by the definition of the triangulated manifold, 3 and 4 
subsimplices are always legal, whilst at this value of /to = there are relatively few legal 
nodes and links (i.e 5-fold coordinated vertices and 4-fold coordinated links). This is to 
be contrasted with fig. f| where the same quantity is plotted for k,q = 2.4. The number of 
legal nodes has increased by nearly a factor of five. 
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In fig. [5] we show the fraction of legal subsimplices P (i) for which an update would 
lead to an geometrically acceptable triangulation. Again the lattices are four dimensional 
with mean volume V = 8000 at zero node coupling. In the case of move (node deletion) 
this is possible with unit probability P (0) = 1, but for subsequent moves P {%) decreases 
until for moves of type 4 (node insertion) it again reaches unity. 

As evidence of the two phase structure mentioned in the introduction we show in fig. 
P a plot of the node susceptibility \ as a function of node coupling kq and for a variety of 
four dimensional lattice volumes V = 500 — 8000. 

X = \ ((iVo 2 ) - (Nof) (12) 

There appears to be a growing peak which shifts and narrows with increasing volume. 
For conventional statistical mechanical systems this would be taken as an indicator of a 
phase transition. This quantity is sensitive to the presence of long range correlations in 
the geometric curvature. The node coupling being (inversely) related to a bare Newton 
gravitational constant, this is taken as evidence that there may be a nontrivial fixed point 
in the theory about which it may be possible to have a consistent quantum theory of 
euclidean gravity. 

Conclusions 

We have described, in some detail, an algorithm to simulate models for quantum gravity 
based on dynamical triangulations. We have demonstrated that the necessary procedures 
can be implemented in such a way that their dependence on dimension is trivial - a single 
compact code can be written in which the dimension is simply input as a parameter. 

Furthermore, we have illustrated the utility of such a code by recording the results of 
some high statistics studies of the four dimensional theory. Numerical evidence is presented 
to support a phase transition. 

SMC would like to acknowledge useful conversations with John Kogut and Ray Renken. 
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Figure 1: Execution profile d = A,V = 8000, k = 



9 



100.0 



80.0 



60.0 



40.0 



20.0 



0.0 



-1.0 0.0 1.0 2.0 3.0 4.0 5.0 

Move type 



Figure 2: Number accepted i-moves per sweep d = 4, V = 8000, k = 
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Figure 3: Number legal subsimplices per sweep d = 4, V = 8000, k = 
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Figure 4: Number legal subsimplices per sweep d = 4, V = 8000, n = 2.4 
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Figure 5: Fraction legal i-moves allowed geometrically, d = 4, V = 8000 
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